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ABSTRACT 

We study diffusion of Cosmic Rays (CRs) in turbulent magnetic fields using test particle simulations. 
Electromagnetic fields are produced in direct numerical MHD simulations of turbulence and used as 
an input for particle tracing, particle feedback on turbulence being ignored. Statistical transport 
coefficients from the test particle runs are compared with earlier analytical predictions. We find 
qualitative correspondence between them in various aspects of CR diffusion. In the incompressible 
case, that we consider in this paper, the dominant scattering mechanism occurs to be the non-resonant 
mirror interactions with the slow-mode perturbations. Perpendicular transport roughly agrees with 
being produced by magnetic field wandering. 
Subject headings: cosmic rays, scattering, MHD, turbulence 
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1. INTRODUCTION 

The interaction between Cosmic Rays, highly energetic 
charged particles, and astrophysical fluids is mediated by 
magnetic fields. As magnetic fields are usually turbulent, 
CRs do not freely stream along these fields but scatter 
(see, e.g., Schlickeiser 2002). Efficient scattering is essen- 
tial for a variety of acceleration mechanisms of CRs, such 
as, for example, diffusive shock acceleration (Krymsky 
1977, Bell 1978, Malkov & Drury 2001 and ref. therein). 

Understanding MHD turbulence is essential for the 
correct description of CR propagation. One popular 
model has been based on the combination of slab and 
two-dimensional perturbations (see Bieber, Smith, & 
Matthaeus 1988). Simplicity of this empirical model 
has appealed to researchers and has been used to ac- 
count for propagation of CRs in solar wind and magneto- 
sphere. Numerical simulations (see Cho & Vishniac 2000, 
Maron & Goldreich 2001, Muller & Biskamp 2000, Cho, 
Lazarian & Vishniac 2002, Cho & Lazarian 2002, 2003), 
however, do not show slab modes, instead, they show 
Alfvenic modes that exhibit scale-dependent anisotropy 
consistent with predictions in Goldreich & Sridhar (1995, 
henceforth GS95) . The scalings of compressible modes is 
still a subject of debate, although it is suggested that 
slow mode is passively advected by Alfven mode (GS95, 
Lithwick & Goldreich 2001), which was verified by nu- 
merics, also fast mode showed relative isotropy which 
was suggestive of a separate acoustic-type cascade (Cho 
& Lazarian 2002, 2003). 

While particular aspects of the GS95 model, e.g. 
the value of the spectral index, has been debated 
(see Boldyrev 2005, 2006, Beresnyak & Lazarian 2006, 
Gogoberidze 2007, Beresnyak & Lazarian 2009a, b), this 
model provide a good start for studying CR scattering. 
This program was realized in a number of publications 
such as Chandran (2000), Yan & Lazarian (2002, 2004, 
henceforth YL02, YL04, respectively), Brunetti & Lazar- 
ian (2007). In the last three papers, following Cho & 
Lazarian (2003), MHD turbulence has been decomposed 
into Alfven, slow and fast modes. 
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In a complex problem of propagation and acceleration 
of CRs we often use so-called diffusive approximation 
which assume that the particle scatter or gain energy 
in small steps. In this approximation the local particle 
dynamics will be averaged to obtain the spatial diffusion 
coefficient, D xx , and the momentum diffusion coefficient, 
-Dpp, that go into the advection-diffusion equation for the 
evolution of quasi-isotropic CR distribution function /: 
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(e.g., Skilling 1975). The source terms has be added 
to the RHS of this equation to account for injection from 
thermal particles and the proper boundary conditions 
should be defined. Here we assumed for simplicity that 
f(x,p) depends only on one spatial coordinate x and the 
magnitude of CR momentum, p. This equation uses "lo- 
cal" system of reference, where particle momentum is 
measured with respect to the rest frame of the fluid. In 
a situation when the advection-diffusion equation is not 
adequate one has to fall back to more general approaches, 
such as Vlasov's equation (see, e.g., Schlickeiser 2002). In 
this paper we study particle dynamics assuming diffusion 
approximation and we monitor if this dynamics looks like 
a diffusive dynamics or not. 

The propagation of CRs is a mature quantitative field, 
which makes use both of analytical studies and numerical 
simulations. For example, a quasi-linear theory (QLT) 
was used to calculate scattering of CRs propagating in a 
mean magnetic field with small perturbations. However, 
as turbulence paradigms were changing, so were the re- 
sults of CR scattering theories. The purpose of this pa- 
per is to measure CR scattering numerically, based on 
the best available direct numerical simulations of MHD 
turbulence and compare these results with what scatter- 
ing theories predict. 

QLT has demonstrated that the gyroresonance in GS95 
type turbulence is substantially suppressed and negligi- 
ble (Chandran 2000, YL02, YL04). However, the key 
assumption of QLT, that the particle's orbit is unper- 
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turbcd, significantly limits its applicability. Addition- 
ally, QLT has problems in treating scattering of parti- 
cles with momentum nearly perpendicular to the mag- 
netic field (see Jones, Birmingham & Kaiser 1973, 1978; 
Volk 1973, 1975; Owens 1974; Goldstein 1976; Felice & 
Kulsrud 2001) and perpendicular transport (see Kota & 
Jokipii 2000, Matthaeus et al. 2003). 

Various non-linear theories have been proposed to im- 
prove the QLT (see Dupree 1966, V61k 1973, 1975, Jones, 
Kaiser & Birmingham 1973, Goldstein 1976). In the re- 
cent paper of Yan & Lazarian (2008, henceforth YL08), a 
nonlinear formalism (NLT) based on V61k (1975) was de- 
veloped. The gyroresonance was found to be marginal in 
incompressible turbulence. However, transit time damp- 
ing (TTD) was fairly efficient which is different from the 
QLT result. TTD due to nonlinear scattering can be 
understood as a scattering by large-scale magnetic com- 
pressions (magnetic bottles formed by slow mode). The 
ideas on perpendicular diffusion has been dominated by 
the field line random walk (Jokipii 1966, Jokipii & Parker 
1969, Forman et al. 1974). It can be justified in a situ- 
ation when CRs do not scatter backwards. However, in 
three-dimensional turbulence, parallel transport is also 
diffusive, and this can reduce perpendicular transport. 

The difference between QLT and NLT has important 
astrophysical consequences. Indeed, in some phases, such 
as hot ISM, the fast mode is strongly damped, which 
leaves only Alfven and slow modes for scattering. Ac- 
cording to the QLT, however, these modes do not pro- 
vide any significant scattering. This would predict that 
there are large volumes in the disk of the Galaxy where 
CRs do not scatter at all, which would question global 
simulations of propagation of CRs in the Galaxy or halo 
made without such an assumption. This will also some- 
what contradict the isotropy of galactic CRs observed on 
Earth, because isotropy suggests efficient scattering. For- 
tunately, NLT corrects this "zero-scattering" QLT pre- 
diction, putting a lower limit on the efficiency of CR scat- 
tering, thus mitigating contradictions described above. 

Test particle simulation has been used to study CR 
scattering and transport before, see, e.g., Giacalone & 
Jokipii (1999), Mace et al (2000), Qin at al. (2002). 
The aforementioned studies, however, used synthetic 
data for turbulent fields, which have several disadvan- 
tages. Creating synthetic turbulence data which has 
scale-dependent anisotropy with respect to the local mag- 
netic field (as observed in Cho & Vishniac 2000 and 
Maron & Goldreich 2001) is difficult and has not been 
realized yet, as far as we know. Also, synthetic data nor- 
mally uses Gaussian statistics and delta-correlated fields, 
which is hardly appropriate for description of strong tur- 
bulence. In contrast, in this paper we are using the re- 
sults of direct numerical MHD simulations as the input 
data for particle scattering simulations. 

Another challenging problem is the back-reaction of 
CRs to the fluid. A particular mechanism for such a 
process, called streaming instability, has been popular in 
describing CR scattering since long time ago (see, e.g., 
Kulsrud & Pearce 1969). In this paper, however, we 
consider only test particle scattering. This describes an 
important physical limit where CR density is negligibly 
small and collective effects are unimportant. Although in 
realistic astrophysical environments CR density is never 
small and CR pressure is always dynamically important, 



the understanding of the test particle limit will give us a 
firm ground for future research into a more general and 
more difficult problem of mutual interaction of CRs and 
MHD fluid. 

In this paper we study the scattering by the incom- 
pressible component of turbulence. If the fast mode is 
present, however, it will dominate scattering of low en- 
ergy CRs, as long as the fast mode is not effectively 
suppressed (YL04). Another very efficient mechanism of 
scattering is the instability between CR fluid and MHD 
fluid (see above). It is present, e.g., when there are strong 
CR gradients that lead to streaming. 

In what follows, we discuss numerical methods, includ- 
ing DNS of MHD turbulence and particle tracing tech- 
nique in § 2. We discuss theoretical expectations for CR 
scattering in § 3. We provide numerical measurements 
of scattering in § 4 and measurements of space diffusion 
in § 5. We discuss our results in § 6. 

2. NUMERICAL METHODS 

In order to trace particle trajectories we were us- 
ing electromagnetic fields obtained in direct three- 
dimensional simulations of MHD turbulence. For the 
purpose of this paper we were using only incompressible 
simulations for a variety of reasons. First, we wanted 
to test those predictions of the theory that pertain to 
incompressible case. Second, the incompressible simula- 
tions were performed with pseudospectral code that has 
explicit dissipation and, unlike finite-difference code has 
no uncertainties due to numerical dissipation. Also, in- 
compressible simulations have larger inertial range. 

2.1. DNS of turbulence 
We solved incompressible MHD equations, 

9 t w ± + S(w T • V)w ± = -z/„(-V 2 )"w ± + f ± , (2) 

written in terms of Elsasser variables which are defined 
in terms of velocity v and magnetic field in velocity units 
b = B/(47rp) 1/2 as w+ = v + b and w = v b, S is 
a solcnoidal projection operator, f ± is Elsasser forcing. 
These arc general equations which can be used for either 
turbulence with no mean magnetic field (i.e. when the 
average of w + — w~ is zero), or in a presence of such a 
mean field. In the latter case perturbations of w + can be 
seen as the waves propagating oppositely to the magnetic 
field direction. Both Alfven and pseudo- Alfven waves 
propagate with the same velocity va = Bq/ {^irp) 1 / 2 . 

We used pseudospectral code described in more de- 
tail in Beresnyak & Lazarian 2009 (a, b) (henceforce 
BL09a,b). The pseudospectral code solves Eq. 2 as 
ordinary differential equation in time for each spacial 
Fourier harmonic, the "pseudo" coming from the fact 
that nonlinear term is calculated in real space, and then 
converted back to Fourier space. The dissipation and 
divergence-free condition for velocity and magnetic field 
are done with simple algebraic operations in Fourier 
space. For time integration we use leapfrog which is 
time-reversible and numerical dissipation is absent, be- 
cause nonlinear term, calculated in this manner, preserve 
both energy and cross-helicity. Therefore the only dissi- 
pation come from the explicit dissipation term. The tur- 
bulence was driven by either independent Elsasser driv- 
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Fig. 1. — Sample spectra from turbulent simulations. Left - 
balanced case, right - imbalanced case. Shown are Elsasser energy 
spectra (see BL09a for more details). 

ing or by pure velocity driving (which formally corre- 
sponds to /+ = / — ). For the purpose of this paper we 
used the results of 768 3 balanced and imbalanced turbu- 
lent simulations from BL09a. Balanced turbulence cor- 
responds to the well-studied limit, where the rms values 
of w + and w~ are equal. Physically this corresponds to 
the situation when the flow of w~ perturbations, which 
propagate along the mean magnetic field direction, bal- 
ances the opposing flow of w + . The more general case of 
imbalanced turbulence is less studied (see BL08, BL09a 
and ref. therein), but more likely to be found in na- 
ture. This is due to the fact, that MHD turbulence is 
often driven by the strong localized source of perturba- 
tions and near the source we mostly see waves moving 
away from the source. Such as a solar wind turbulence 
near the Sun, which is strongly imbalanced 3 . Naturally, 
we are also interested in particle scattering in the imbal- 
anced turbulence, although few theoretical predictions of 
scattering exist in this case, if any. 

Another dimension in parameter study of BL09a was 
the strength of perturbations with respect to the mean 
field. The SB ~ Bq is called trans- Alfvenic case, where 
perturbations are of the order of the mean field. We 
also consider sub- Alfvenic case when perturbations were 
approximately 10 times weaker than mean field Bq, 
which correspond to so called Alfvenic Mach number 
Ma ~ 0.1. The latter case can be considered as smaller 
scales of trans- or super-Alfvenic turbulence that can- 
not be reached directly by 3D simulations of aforemen- 
tioned flows. In order for sub- Alfvenic turbulence to be 
strong 4 it has to be driven anisotropically on its outer 
scale, which was realized in BL09a,b. The computational 
box was also elongated in the direction of the mean field, 
with parallel size 10 times larger than perpendicular size 
for Ma — 0.1 case. Throughout the paper, when we men- 
tion the "box size" and the "outer scale of turbulence" it 
means perpendicular size, the parallel size is the same for 
trans- Alfvenic cubes and 10 times larger for sub- Alfvenic 

3 The measurement of the imbalance in the solar wind has been 
possible with the advent of satellites that independently measure 
the velocity and the magnetic field at the same point. Similar 
measurements for other astrophysical sources, such as ISM, are yet 
to be developed. 

4 Strong MHD turbulence appears naturally as a result of the 
anisotropic cascade. Even if turbulence is driven weakly with re- 
spect to the mean field, the perpendicular cascade of weak turbu- 
lence (Galtier et al, 2000) will increase the strength of interaction 
until it becomes strong. The realistic ISM turbulence, however, is 
driven strongly, such as SB ~ Bq on the outer scale. So turbulence 
is strong to begin with and continue to be strong along the cascade 
(GS95). 



cubes. 

The note of caution has to be said with regard to 
sub- Alfvenic simulations being the small scales of trans- 
Alfvenic turbulence. As we use periodic boundaries, the 
scales which are larger than the cube size are excluded 
from consideration. This means that we cut out a range 
of scales in the inertial interval of turbulence and all 
larger scales are represented only by the value of the 
mean magnetic field (the mean velocity can be excluded 
by the local frame of reference). This could or could 
not be satisfactory for simulations of particle scattering. 
If the resonant scattering mechanism is effective, then 
particles mostly interact with those scales of magnetic 
perturbations that are present in the simulation. In the 
opposite case the aforementioned interaction can be less 
effective than the interaction with large scale perturba- 
tions that are not present in the numerical cubes. In this 
case the result should not be trusted. We will return 
to this question below. Accidentally, the QLT consider 
scattering in a manner which is consistent with approach 
that ignores larger scales and consider particle gyrating 
along a strong guiding field and interacting with small 
resonant perturbations. 

The turbulence was driven with specially designed 
quasi-stochastic driving on outer scale (with wave num- 
bers in the interval k = 2. .3. 5) with self-correlation time 
of 2 in code units. The driving worked until stationary 
state was reached. The scale of the largest coherent eddy 
in the simulation was around 0.2 of the cube size. This 
is the outer scale of turbulence L = 0.2. This largest cor- 
relation scale of velocity and magnetic perturbations is 
determined by nonlinear interaction and is typically less 
than the driving correlation scale. On the driving scales, 
i.e. k — 2. .3. 5 the turbulence is not yet fully developed 
and the spectrum is distorted, having a characteristic 
bump, and well-developed turbulence starts with k = 5. 
Another definition of outer scale is through anisotropy. 
One can expect the anisotropy to follow a Goldreich- 

Sridhar critical balance fcii = k 2 / 3 L^ 1 / 3 . This also lead 
to the estimate of L = 0.2[] 5 , with turbulence being ap- 
proximately isotropic at k — 5 and approximately 1 : 2 
anisotropic at k = 40 for trans- Alfvenic case. At any 
given time the cube contained large number of indepen- 
dent turbulent realizations (> 40). Spectra for one bal- 
anced and one imbalanced case are presented on Fig. [TJ 
Further details of the code and simulations can be found 
in BL09(a,b). 



2.2. Particle tracing 

The electric field in the laboratory frame was obtained 
through E = ~[vx B]/c equation, assuming Va/c = 10~ 5 
(a typical value for the ISM) 6 . The particles were injected 
randomly through the cube and the trajectories were 
traced by hybrid Runge-Kutta quality-controlled ODE 
solver, assuming periodic boundaries for particles as well 
as fields. 

In particular, we solved 6 equations: 

5 Here we omit 2n factor normally present in size-wavevector 
relation, as we normalized cube size to unity. 

6 The velocity was measured in the Alfvenic units and the electric 
field is in the same units as magnetic field 



Fig. 2. — Diffusive behavior of the particles in the tracing simu- 
lations. For different we show ensemble- averaged square devi- 
ations, which are proportional to time, x and y are measured in 
units of the cube size. 



Fig. 3. — Cartoon illustrating that the CR gyroresonance scat- 
tering in strongly anisotropic (GS95) turbulence is very ineffective. 
In perpendicular direction CRs feel many uncorrelated eddies and 
the interaction averages out. 
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Here u is the normalized space component of the 4- 
velocity, u = u/70, where 70 is the initial particle 
gamma- factor. Also 7 = \J 1 J~ffi + u 2 . s — (eB /mc 2 )s 
is a self-time measured in cyclotron frequency units (a 
gyration frequency in particle's own frame). A particle 
with /i = will make a full orbit in Bo field in 2ir time. 
Therefore, we conveniently measure scattering frequency 
relative to gyration frequency. The measure of initial 
particle's energy, normalized Larmor radius is expressed 
as vl = mc 2 7o/e_Bo. Physically, one can think of 70 is a 
measure of relativicity of the particle, i.e. for small 70 we 
will recover nonrelativistic equations, and for large 70 - 
ultra-relativistic equations. At the same time, r*/, is the 
measure of energy, but with respect to the perpendicu- 
lar size of the simulation box. In most simulations we 
took I/70 (which enters only in the equation for 7) as 
zero or close to zero, such as 10~ 5 , this corresponds to 
ultra-relativistic particles. The rx was varied from 0.1 
of the cube size to around a grid size. Fig. [2] presents 
ensemble-averaged square distance vs time for different 
ri,. The square distance grows linearly with time, which 
is expected for diffusive motion. 

3. EXPECTED CR TRANSPORT PROPERTIES IN MHD 
TURBULENCE 

3.1. Formalism for NLT 

We start with explaining QLT which is the theory 
for resonant interactions: gyroresonance scattering and 
transit scattering (also called transit time damping, 
TTD). The resonant condition is lo — kuvfj, = nVt (n = 
0, ±1,2...), where u is the wave frequency, SI = fio/7 is 
the relativistic gyration frequency, fi — cos 0, 8 is the 
pitch angle of particles. TTD corresponds to n = and 
it requires compressible perturbations. Most of the gy- 
roresonance contribution comes from n = 1. 

It was demonstrated that scattering by Alfvenic tur- 
bulence is substantially suppressed due to its anisotropy 
(Chandran 2000, YL02). Fig. H illustrates why interac- 
tion is suppressed. The scattering rate in GS95 turbu- 
lence with outer scale of L and assuming that 9 is not 
close to is given by QLT as (YL02): 
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where T[a, z] is the incomplete gamma function, k max 
correspond to the dissipation scale of turbulence and 
f t \\,res v f 1 — ^- The scattering frequency, therefore, is ap- 
proximately Bohm-likc if Larmor radius is of the order of 
L, but then falls steeply as f2~ 15 and becomes negligible 
for small energies. 

Contrary to QLT which assume that the magnitude 
of the magnetic field stay constant, NLT relaxes this as- 
sumption and allow this quantity to change gradually, 
adiabatically with respect to particle motion. Due to 
conservation of adiabatic invariant p\ / B (see Landau 
& Lifshits 1975) the pitch angle will gradually vary, re- 
sulting in resonance broadening (Volk, 1975). Nonlin- 
ear transport (NLT) formalism is based on the replace- 
ment of the sharp resonance between waves and particles 
<5(fc||Z)|| — uj ± nfi) from QLT to the "resonance function" 
R n (YL08): 



Rn, = U / dte 



£(fc|t zj|| +rof2— uj)t— ^ fcii Vj_t 



|fe||Au|| 



■ exp 



(knvfi — cu + nCiy 



(6) 



The width of the resonance function depends on the 
perturbation strength of the turbulence A/i = Avu/vj_ ~ 

\J5B/B = \JMa)- For gyroresonance (n = ±1,2,...) 
the result depends on whether fi is strongly or weakly 
perturbed by regular field. If fi » A^, the result is 
similar to QLT, because the exponents in Eq.© become 
close to ^-functions. For /i < A//, however, the result is 
different. To demonstrate this we can consider the case of 
90° scattering. Indeed, if /1 — > 0, the resonance happens 
mostly at k\\ res ~ Q/Av, while in QLT fcii res ~ f2/uii — > 
00. Compared to TTD, however, for gyroresonance 
is still smaller in incompressible case (see Fig. 0|) due to 
anisotropy. 

3.2. Scattering in strong MHD turbulence 

Assuming the tensor of magnetic perturbations intro- 
duced in Cho et al. (2002), which is consistent with the 
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Scattering in strong incompressible trubulence 



10" 

o 

cc 10 

a 

10" 12 



10" 




— TTD 
- - gyro Alfven 
--gyro slow 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Fig. 4. — Scattering of CRs with R = r^/L = 10 in strong in- 
compressible trans- Alfvenic turbulence. Solid line represents TTD 
contribution. Dash-dot refers to the gyroresonance with the Alfven 
mode while dash line is for the gyroresonance with the slow mode. 
See eqs [7] also see YL08. 

GS95 model, we can calculate scattering from TTD and 
gyroresonance. Assuming small energies corresponding 
to Larmor radii much smaller than the outer scale one 
gets for TTD (YL08): 
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where £ = kL, Ei(£) = dtexp(—£t)/t and q = 

{U,ma X Ml)-*/\ 

Fig. 2] displays the pitch angle diffusion coefficients re- 
sulting from TTD scattering and gyroresonance (YL08). 
We see that gyroresonance is mostly subdominant, how- 
ever at small pitch angles TTD is inefficient and gyrores- 
onance dominates. 

4. MEASUREMENTS OF SCATTERING 

-D/i/i scattering property was measured in the trac- 
ing experiments where an ensemble of particles with the 
same Tl (energy) and a particular /j,q were traced by a 
certain time. This time was determined by the condition 
that the rms of deviations of /i is small (i.e. 0.1-0.01). 
Then the curves of the ensemble- averaged ((/z — Mo) 2 ) 
were fitted with a linear curve, and so D MM was obtained. 
The for trans- Alfvenic case is presented on Fig. [5l 
For sub-Alfvenic case we noticed that there were very 
few 90° scattering events. This will be explained in the 
next section. 

As we see from Fig. [5] the measurement of scattering 
frequency is incompatible with QLT. The scattering fre- 
quency normalized to the gyration frequency is propor- 
tional to the Larmor radius i.e. it is constant with en- 
ergy (as Ctri, = v « c). It would be reasonable to assume 
then that particles of all energies scatter on the same ob- 
jects, magnetic bottles, formed by large scale slow-mode 
perturbations. The same result could be obtained from 
NLT, taking into account A/i ~ \x in strong turbulence. 
At larger energies scattering becomes less efficient i.e. 
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Fig. 5. — Dimensionless CR scattering coefficient D MM /S7 for the 
case of fi = 0.71 vs Larmor radius rt expressed in cube size units 
(solid line). We suppose that it is dominated by magnetic bottles 
formed by slow mode, this is why D^^ (dimensional scattering fre- 
quency) is almost constant, i.e. independent on particle's energy. 
For comparison, we plot various theoretical predictions: QLT pre- 
diction for Alfven and slow mode (dashed); QLT prediction for fast 
mode (dot-dashed), note that in our data fast mode was absent, 
so this line is only for reference; hypothetical Bohm scattering or 
maximally efficient scattering (dotted). 

high energy particles "feel" less mirrors. This transition 
happen at around Tl/L ~ 0.1. 

Due to the lack of outer scale, the Z? MA1 for sub-Alfvenic 
case is supposed to be QLT-like and very small. As such 
it was severely contaminated by numerical error, in par- 
ticular, fields interpolation error and was not obtained in 
this study. 

5. MEASUREMENTS OF SPACE DIFFUSION 

The measurements of space diffusion D xx and D yy were 
more straightforward than the measurements of be- 
cause we did not limit the integration time as in the 
previous section. Therefore, we integrated for as long 
as it took for the ((x — xo) 2 ) and ((y — yo) 2 ) to show 
good diffusive linear dependence with time. Those inte- 
gration times turned out to be very long, so the particles 
crossed outer-scale of turbulence for many times. There- 
fore these measurements correspond to diffusion on outer 
scale and not to "sub-diffusion" (see, e.g., YL08). More- 
over, we measured diffusion with respect to some global 
frame of reference, determined by the global mean mag- 
netic field. Therefore, our measurements do not neces- 
sarily correspond to theories that measure "parallel" or 
"perpendicular" diffusion with respect to the magnetic 
field lines. 

Also, we were only able to obtain the lower limit of 
D xx for sub-Alfvenic case due to the very low 90° scat- 
tering frequency. This was manifested by the fact that as 
we increased the precision of the code, the D xx increased 
and did not show convergence. This very low scattering 
frequency has to do with what we discussed earlier - the 
lack of larger scale perturbations in sub-Alfvenic cubes. 
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Fig. 6. — Parallel diffusion in trans-Alfvcnic case. Solid — "bal- 
anced" , dotted - " slightly imbalanced" dashed - " strongly imbal- 
anced", dot-dashed — "very strongly imbalanced" (see BL09a for 
more details). 

In this case, since A/x < /i, the resonance function be- 
comes narrow so that marginal interaction is available at 
90° 7 . In other words, in the absence of the large-scale 
perturbations, which are normally present in nature, but 
absent in our sub-Alfvenic cubes, the 90° scattering be- 
comes problematic and parallel diffusion is replaced by 
ballistic propagation along mean field. At the same time, 
this suggests that QLT (resonant scattering) can not be 
used for low energy particle scattering, as large scales 
contribute more than the resonant scales. 

As the mean magnetic field was along 'x' axis, our D xx 
coefficient correspond to "parallel diffusion" , while D yy 
correspond to "perpendicular diffusion" . This correspon- 
dence, however, is tentative, since most theories predict 
"parallel" or "perpendicular" diffusion as happening with 
respect to the local magnetic field lines. We nevertheless 
will use terms "parallel" and "perpendicular" to D xx and 
Dyy. We also claim that the measurement of the diffu- 
sion with respect to the global reference frame has more 
practical importance and is easier applicable to the re- 
sults of observations. 

5.1. Parallel diffusion 

The results for parallel diffusion for trans- Alfvenic case 
are presented on Fig. [6] Along with standard 'bal- 
anced' MHD turbulence case (presented by solid line) we 
calculated this diffusion coefficient for simulations with 
different degree of imbalance, using datacubes from sim- 
ulations of Beresnyak & Lazarian 2009a. As the afore- 
mentioned paper (along with the earlier study Beresnyak 
& Lazarian 2008) are the first high-resolution simulations 
of stationary strong imbalanced turbulence, it is impor- 
tant 8 to numerically study the scattering coefficient, even 
more so when the theory is lacking. 

As we see from Fig. [5] at small energies D xx /Q, is lin- 
early proportional to ri,, i.e. as in the case with D^ 
the scattering frequency is independent of energy. At 
higher energies it becomes proportional to the square of 

7 The case with sub-Alfvenic turbulence may be viewed as a 
transition to the QLT case, where the resonance function shrinks 
to S function and there is no TTD resonance at 90° 

8 For example, Solar Wind exhibit imbalanced turbulence up to 
distances of around 1 AU, with perturbations coming from the Sun 
being prevalent over backward-going perturbations. 
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Fig. 7. — Perpendicular diffusion in sub-Alfvenic case SB/B = 
1/10 

Tjj. This is again consistent with the behavior of D^^/fl 
from Fig. i.e. that at high energies ~ fl and as 
D xx ~ 1/Dpfj,, D xx /Vt ~ 1 / ri 2 ~ r\. The transition hap- 
pens at Tl/L 0.1, same as in Fig. [5j So we conclude 
that the measurements of D xx and are consistent 
with each other. A note of caution towards direct com- 
parison of these two measurements is due, however. In 
the measurement of D xx we did not control particle's 
energy, which could undergo changes during the long in- 
tegration times of D xx measurement. D^, however, was 
measured during short times, and as electric field was as- 
sumed small (smaller than B by a factor of va/c ~ 10~ 5 ), 
there wasn't significant energy change during this short 
time. This can explain why the transition between two 
regimes of scattering is more sharp of Fig. [5] rather than 
Fig. [6l Also, as we mentioned previously, D xx is the dif- 
fusion coefficient measured in the global reference frame, 
while -D MM defines pitch-angle scattering with respect to 
the local field direction. 

Fig. [5] also shows that the diffusion coefficient is pretty 
independent on the degree of imbalance, indicating that 
the trans-Alfvenic imbalanced turbulence has approxi- 
mately as many magnetic bottles as its balanced coun- 
terpart. This is consistent with the assumption that only 
large-scale perturbations significantly contribute to scat- 
tering. Indeed, in the imbalanced simulations of BL09a 
the outer-scale magnetic field was determined primar- 
ily by the stronger Elsasser component and has similar 
structure and magnitude and outer-scale magnetic field 
in balanced simulations. 

5.2. Perpendicular diffusion 

Perpendicular diffusion coefficients are presented on 
Figs. [3 [SI As to various models, regimes and termi- 
nology of perpendicular diffusion we refer the reader to 
YL08 and refs therein. Let us first interpret the mea- 
surements in the sub-Alfvenic case. We chose initial par- 
ticle's pitch angle to be 45°. As we discussed earlier, due 
to the particular choice of the data, there wasn't any 90° 
scattering in this case. I.e. particles moved ballistically 
along 'x' axis, but their trajectories diffused from the 
center due to magnetic field wandering. As suggested by 
Fig. [7] the dependence of this plot is almost linear, i.e. 
Dyy /CI r/L or D yy is independent of energy. 

In the three-dimensional turbulence, field lines are 
diverging away due to shearing by Alfven modes (see 
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Fig. 8. — Perpendicular diffusion in trans-Alfvenic case SB/B = 

1. 

Lazarian & Vishniac 1999, Narayan & Medvedev 2002). 
Most recently the diffusion in magnetic fields was con- 
sidered for thermal particles in Lazarian (2006, 2007). 
The cross-field transport can result from the deviations 
of field lines at small scales, as well as field line random 
walk at large scale. 

If we assume that the particle follow magnetic field 
line and is diffused only by the outer-scale magnetic field 
wandering, the perpendicular diffusion can be expressed 
as D yy /(L 2 n) a 2- 1 / 2 (L tr / L) 2 ■ (L/L\\) ■ (r L /L), where 
JDii is the outer parallel scale (which is 10 times bigger 
than L in our sub- Alfvenic simulation) , 1 / v2 is the co- 
sine of pitch angle and L tr is the distance the particle is 
deflected when it travels Ln along the field line (see eq. 
(26) in YL08). We would expect L tr to be close to L. 
From the fit of Fig. [7] we derive L tr /L ss 0.92 which is 
fairly close, considering the uncertainty in L. Using the 
same argumentation we obtain L tr /L ss 0.53 from the 
fit of Fig. |6] which is short of what we expected. This 
is an indication that the impediment of travel in parallel 
direction which is present in trans-Alfvenic case due to 
90° scattering decreases diffusion in perpendicular direc- 
tion. We stop with this conclusion, as there is clearly 
not enough data for a detailed comparison with different 
models in YL08. 

6. DISCUSSION 

In this paper we numerically measured diffusion coeffi- 
cients that arise when particle propagates in a turbulent 
magnetic fields. Unlike previous studies, we used realis- 
tic fields obtained in a three-dimensional simulations of 
MHD turbulence. The focus of this paper was the incom- 
pressible case, where fast magnetosonic mode is absent. 
The earlier QLT calculations presumed that particle scat- 
tering is negligible in this the perturbations are 
extremely anisotropic with respect to the mean field. We 
figured that QLT is not applicable when the magnitude 
of the magnetic field is strongly perturbed, and that an- 
other approach called NLT has to be adopted (YL08). 
NLT allows relatively efficient scattering through TTD 
as the particle's pitch angle changes adiabatically and 
makes possible for 90° scattering. One can interpret this 
as scattering through large-scale magnetic mirrors. 

We confirmed this picture of mirrors by measuring 
scattering frequency which is independent on energy, for 
small energies. We also studied spacial diffusion which, 



in the case of parallel diffusion is related to scattering 
frequency. The case of perpendicular diffusion is more 
complicated. We showed that if particles do not scatter 
in parallel direction, the perpendicular diffusion is mostly 
due to magnetic field line wandering. This case could be 
unphysical though, as the absence of parallel scattering 
was due to the absence of larger scales in the simulation. 
In the case when parallel diffusion was operating, the 
perpendicular diffusion was reduced. At this point, how- 
ever, we don't have enough data to distinguish between 
different models of perpendicular diffusion. 

The special attention should be brought to the astro- 
physical interpretation of scattering in the imbalanced 
turbulence. Figs. 6-8 indicate that the scattering is simi- 
lar to the balanced case. This is qualitatively and quanti- 
tatively agree with the picture that was presented in this 
paper, namely that in the incompressible case most of the 
scattering will come from the outer scale of turbulence 
and most of the perpendicular diffusion will come from 
the field wandering on the outer scale. This fact, how- 
ever, does not mean that the scattering in the astrophys- 
ical objects will be the same regardless of the degree of 
imbalance. The key to understand this is to understand 
the nature of our MHD simulations. In these simulations 
we kept the fluctuation amplitude and the anisotropy 
controlled on outer scale, the physically all-important 
dissipation rate, however, varied greatly depending on 
the degree of imbalance. In astrophysics, turbulence is 
caused by the sources of kinetic energy, such as stellar 
and AGN jets, stellar winds interacting with the ISM, 
supernovae, the Sun creating perturbations in the solar 
wind. Turbulent dissipation will have to balance this in- 
flux of energy, however, dissipation depends greatly on 
the degree of imbalance, therefore, in a situation with a 
constant influx of energy imbalanced turbulence will have 
much larger perturbation amplitudes, which will result in 
a much more efficient scattering. With respect to relation 
between dissipation rate and perturbation amplitude we 
refer to the imbalanced turbulence model presented in 
Beresnyak & Lazarian 2008, as the most realistic model 
to-date, and the simulations in BL09a. A word of cau- 
tion towards directly using these results is due, however. 
Aforementioned model and the simulations in BL09a de- 
scribe stationary imbalanced turbulence. However, as 
we learned from these studies, the time of establishment 
of stationary state greatly increases with larger imbal- 
ances. As astrophysical processes are usually transient, 
it is possible that in a situation with large imbalance 
the stationary state will not be achieved. The station- 
ary imbalanced turbulence could still be used to infer the 
propertied of small-scale fluctuations, as timescales are 
smaller, but, as we saw in this study, large scales are im- 
portant for scattering. This problem will be solved by 
the models of transient and inhomogeneous imbalanced 
turbulence, although at present such models are still in 
their infancy (see, e.g., the Appendix in BL09a). We are 
optimistic, however, that the properties of CR scattering 
in realistic astrophysical objects that feature imbalanced 
turbulence, such as solar and stellar winds, AGN jets and 
many others will be figured out. 

Although NLT prediction for nonresonant mirror scat- 
tering by incompressible turbulence puts a lower limit 
on CR scattering, in most realistic astrophysical circum- 
stances, there several mechanisms that could compete 
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with it. If the fast mode is present and have a sufficient 
amplitude in the range of scales corresponding to Lar- 
mor radii of low-energy CRs, those CRs will be scattered 
primarily due to fast mode (YL04). Also, if CRs have 
large density gradients and tend to stream in a particu- 
lar direction, such as CRs escaping the Galaxy, or CRs 
streaming in front of the supcrnovae shock, the back- 
reaction to MHD fluid will be important. Speaking of 
turbulence, in this paper we considered generic astro- 
physical turbulence driven on large scales. Other types 
of astrophysical turbulence are often important for scat- 
tering and acceleration. This includes MRI turbulence 
(see, e.g., Hawley et al, 2001), and turbulence generated 
by CR-MHD fluid interaction in supernova shocks (see, 
e.g., Beresnyak et al 2009 and rcf. therein). 



Stochastic acceleration by MHD turbulence was not 
studied here as the correct calculation requires time- 
dependent MHD fields, so that simulations of turbulence 
are integrated at the same time as particles propagate. 
This will be a matter of a future research. 
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